Nematic phase without Heisenberg physics in FeAs planes 
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We use Monte Carlo simulations and analytical arguments to analyze a frustrated Ising model 
with nearest neighbour antiferromagnetic coupling Ji and next nearest neighbour coupling J2. The 
model is inspired on the physics of pnictide superconductors and to some extent we argue that it 
can be more representative of this physics than the Heisenberg counterpart. Parameters are chosen 
such that the ground state is a columnar or striped state, as observed experimentally, but is close 
to the transition to the simple Neel ordered antiferromagnetic state R = J2/IJ1I > 0.5. We find 
that a nematic phase is induced close to 7? = 0.5 by finite size effects and argue that this explains 
experiments in imperfect samples which find a more robust nematic state as the quality of the 
sample decreases [A. Jesche et al, Phys. Rev. B 81, 134525 (2010)]. Including the effect of a 
weak coupling with the lattice we find that a structural transition occurs associated with a nematic 
phase, with a magnetic transition occurring at a lower temperature. These two transitions merge 
into a single structural and magnetic transition with a stronger first-order character for larger spin- 
lattice couplings. These two situations are in agreement with the different phenomenologies found 
in different families of pnictides. 

PACS numbers: 74.70.Xa; 75.10.-b; 75.40.Mg; 75.40.Cx 



I. INTRODUCTION 

The FeAs high-Tc materialsii^ are a very interesting 
new playground to study the interplay between lattice, 
magnetism and superconductivity. Quite generically un- 
doped or slightly doped samples show a magnetic phase 
which breaks the C4 symmetry of the lattice, consisting of 
magnetic moments aligned ferromagnetically on one di- 
rection and antiferromagnetically in the other direction. 
This "columnar" or "spin-stripe" phase is accompanied 
by an orthorhombic distortion of the high temperature 
tetragonal lattice. 

From symmetry considerations it is quite natural to 
assume that the structural transition from the high tem- 
perature tetragonal phase to the orthorhombic phase is 
driven by the magnetism. However '1111' materials like 
LaOFeAs'^ display the structural transition at a higher 
temperature than the magnetic transition. This has led 
some authors to propose that the structural transition 
drives the magnetism. Even more surprising is the fact 
that the temperature splitting between structural and 
magnetic transitions decreases upon increasing the sam- 
ple quality^ and the transport anisotropy is reducedi^ 

The stripe-magnetic phase, breaking both the C4 sym- 
metry and the translation symmetry can be called "smec- 
tic", in analogy with the crystal phase that elongated 
molecules form in liquid crystals. However in the last 
few years the question has emerged whether a "nematic" 
phase can also occur in these systems. The term "ne- 
matic", which is also borrowed from the field of liquid 
crystals, indicates a phase where the rotational symme- 
try is broken, while the translational symmetry is fully 
preserved. In the pnictide case, for instance, the square 



lattice C4 symmetry of the FeAs layers in the tetrag- 
onal phase could be reduced to C2 by the occurrence 
of a purely magnetic nematic phase. In turn this ne- 
matic state could involve an interplay between struc- 
tural and magnetic properties giving rise to structural 
and lattice signatures even in the absence of a "smectic" 
magnetic order. This issue not only has been addressed 
from the theoretical point of view^--, but increasing evi- 
dence is now experimentally emerging of a nematic state 
above the magnetic transition in pnictides. Scarming tun- 
nelling microscopy studies detect quasiparticle electronic 
states having a reduced C2 symmetryi^i as well as lo- 
cal anisotropies^^, transport experiments in detwinned 
122 crystals find an in-plane anisotropy, which cannot 
be attributed to lattice distortions onl y^^'^^ , while opti- 
cal experiments find anisotropic mid-infrared spectrapi^. 
If such nematic phase exists, the coupling between the 
magnetic nematic state and the lattice would naturally 
account for the occurrence in some materials (like 1111 or 
Ba(Fei_a;Co2:)2As2) of the tetragonal-orthorhombic tran- 
sition at a higher temperature than the magnetic tran- 
sition. The coupling between nematicity and lattice has 
also been taken as a possible explanation for changes in 
elastic properties observed by ultrasound spectroscopy in 
FeAs systems-^ , and the suppression of orthorombicity in 
the superconducting phase^^. 

An important question is if an electronic nematic phase 
drives the lattice distortion or vice versa. The split tran- 
sitions occur only when the thermal transition is con- 
tinuous or very close to contirmous. On the other hand 
most '122' materials like SrFe2As2i^ show a first-order 
transition in which simultaneously the structural and the 
magnetic order parameter become non-zero. 
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The more straightforward explanation for nematic 
phases starts from a two-dimensional frustrated Heisen- 
berg model^ and is based on theoretical results that dates 
back two decades.-^ This explanation relays on the order 
by disorder mechanism by which the degeneracy of the 
frustrated Heisenberg model at the classical level is bro- 
ken by thermal and quantum fluctuations. As a conse- 
quence a nematic phase is stabilized which has Ising sym- 
metry. In two dimensions antiferromagnetic (AF) long- 
range order can only occur at zero temperature while 
the nematic phase, belonging to the Ising universality 
class, can occur at finite temperatures thus a large re- 
gion of temperatures exists where the system is in a ne- 
matic phase. Interlayer coupling stabilizes a three dimen- 
sional AF state but always with a magnetic transition 
occurring at a lower temperature than the Ising nematic 
transition!^ 

While this explanation is very appealing it is worth 
to examine whether the Heisenberg physics is really nec- 
essary and if there can be other mechanism by which 
nematic phases can be stabilized. One motivation to do 
so is that the region where the Heisenberg physics is rel- 
evant can be quite narrow around the Neel temperature. 
Indeed stoichiometric FeAs materials have small mag- 
netic moments and are metallic which indicates rather 
itinerant magnetisnti^. The magnetism thus is rather 
collective which means that a given site does not inter- 
act only with the neighbors but with several sites over 
a coherence distance 3> a, with a the lattice spacing. 
To observe Heisenberg-like critical fluctuations the mag- 
netic correlation length ^ has to exceed which occurs 
only very close to the magnetic critical temperature. If 
^0 is sufficiently large a small Ising like anisotropy or a 
three dimensional coupling will make the 2D Heisenberg 
physics and the Goldstone modes rather irrelevant. One 
way to suppress from the start the presence of the Gold- 
stone modes and have a finite critical temperature is to 
consider Ising spins. 

In this work we study a frustrated two-dimensional 
Ising model. For simplicity we assume localized 
spins which interact through effective nearest-neighbor 
and next-nearest-neighbor antiferromagnetic interactions 
(Ji > and J2 > respectively) (see, e.g. Ref. We 
thus obtain a frustrated Ising model without any inter- 
action between the spins and the lattice and only at a 
later stage we will introduce a spin-lattice coupling. The 
model's phase diagram is studied using Monte Carlo sim- 
ulations and finite size analysis. 

This model is interesting in its own and although 
its magnetic properties have already been extensively 
investigated^i, our analysis will evidentiate new features 
and it will find interesting connections with the real pnic- 
tide systems. We show that the critical behavior corre- 
sponds closely to a 4-component Potts model. 

For < R = J2/\Ji\ < 0.5, it is well known that 
there is a second-order transition toward a low temper- 
ature Neel antiferromagnetic phase. For R > 0.5 the 
system shows a magnetic transition toward the "spin- 



stripe" phase observed in FeAs planes and mentioned 
above. While in other cases, like in manganites, this 
antiferromagnetic spin configuration has been named c- 
type antiferromagnetism, we will keep the "spin-stripe" 
phase nomenclature as usual in the pnictide field. 

Our analysis aims to explore the possibility of hav- 
ing nematic phases when all the Heisenberg physics is 
suppressed. Although the model is very simplified we 
show below that it explains well several aspect of the 
phenomenology of FeAs planes. In particular we find 
that nematic phases can be stabilized at finite tempera- 
tures by finite size effects. We argue that such effect can 
explain the larger tendency of "bad samples" to display 
nematic phases. 

The paper is structured as follows. In Sect. II we 
introduce the model and the technical framework to solve 
it. In Sect. HI we report the results in the absence of 
the spin-lattice coupling, which is instead introduced in 
Sect. IV. Our concluding remarks are in Sect. V. 

II. MODEL AND METHODS 

We consider a two-dimensional frustrated Ising model 
with antiferromagnetic nearest-neighbour (nn) and next- 
nearest-neighbor (nnn) interactions. The Hamiltonian of 
the model reads 

N 

i—1 S—x.y 
N 

+ 51 J2crzcri+v, (1) 

with Ji and J2 both positive. The notation i-\-x {i+x+y) 
indicates the first (second) neighbour of site i in the x 
{x -f y) direction. One can use a canonical transforma- 
tion (Ti — >■ {—ly^'^y^ai to change the sign of the first 
term. This means that since the model is classical it is 
fully equivalent to a nn ferromagnetic model with AF 
frustration among nnn which has been extensively stud- 
ied in the literature'^^. 

If frustration is strong enough {R > 0.5), the ground 
state is given by a spin-stripe configuration. This configu- 
ration breaks both rotational symmetry and traslational 
symmetry in one direction. This state is characterized 
by a staggered magnetization with (tt, 0) or (0, tt) wave 
vectors, which, for one configuration, is defined by 

"^(o.-) = ^E^^(-l)'^ (2b) 

j 

We notice that at T = 0, for a given orientation of the 
spin stripes, only one of the two staggered magnetizations 
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is non zero. The sign of the two parameters depends on 
the position of the spin up stripes (odd or even rows 
or columns). In Fig [T] we report the phase diagram of 
the model in Eq.(IT]) as obtained by Moran-Lopez et al^ 
studying the equivalent nn ferromagnetic version. Here 
we have added two points obtained from our MC calcula- 
tions. More detailed data including nematic and lattice 
effects will be shown below. In Ref. 21, it was found that 
the transition toward the spin-stripe phase is weakly first- 
order up to a tricritical point at i? = 1.144. For higher 
values of R the transition becomes of the second-order. 




Stripe-AF 



Figure 1: (Color online) Phase diagram of the system. Each 
point represents the critical temperature for a given value of 
R, and the types of point depends on the different meth- 
ods used to compute them. Empty circles (green online) 
have been obtained with Monte Carlo; empty triangles (blue 
online) have been obtained with real space renormalization 
group; crosses (purpleonline) indicate have been obtained 
with Monte Carlo renormalization group; the solid (yellow 
online) circles have been obtained with series expansion. The 
open square (light blue online) marks the tricritical point at 
R = 1.14421. The solid square (red online) mark the critical 
temperature values obtained with our MC simulations. The 
model studied in Ref. [U considers ferromagnetic nn interac- 
tions so that the ordered phase up to _R = 0.5 should more 
precisely be named "Neel-F". 

In Fig. [5] we show a domain of the perfect stripe order 
in which the directionality is preserved but the phase of 
the AF changes by tt. At finite temperatures it can hap- 
pen that domains of this kind are preferentially excited 
destroying the long-range magnetic order but still break- 
ing C4 rotational symmetry so that the system reaches a 
nematic state. 

Following Ref. d we define a nematic magnetization 
which measures the global breaking of C4 symmetry and 
which for one configuration reads, 



N 



Here the local nematic parameter 
01 = 7:0-i{<^i+'. 



is defined by 



H+y) 



(3) 



(4) 



The magnetic and nematic order parameters are defined 
as the thermal averages (...) of Eqs. (gal , ([23) ,([31) . The 
range of all these three parameters is [—1 : 1]. 
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Figure 2: (Color online) Illustration of a spin domain (gray 
arrows, red online), in which the staggered magnetization in 
X direction has opposite sign with respect to the embedding 
cluster (black arrows). 
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The "nematic phase" is characterized by finite values of 
the nematic magnetization, consequence of the break of 
the rotational symmetry C4, but vanishing values of the 
staggered magnetization. It is customary to call (m^) 
the nematic order parameter— and we follow this conven- 
tion here, although it is different from zero also in the 
spin-stripe or smectic phase which, strictly speaking, are 
not nematic phases. The sign of (m^) depends on the 
preferred orientation of the spin stripes. Table |lT] sum- 
marizes the above discussion of the order parameters. 

To numerically implement our analysis we define a 
square matrix (that is our lattice) , every element of which 
may assume the values ±1, according to the up or down 
orientation of the spin on the site. Physical quantities 
like staggered magnetizations, nematic order parameter, 
spatial correlation functions and susceptibilities, are cal- 
culated using the MC method^^-^'^. This method exploits 
the Metropolis algorithm to generate, in a random way, 
a chain of states (called "Markov chain"), which are dis- 
tributed according to the Boltzmann distribution func- 
tion. Thus the MC procedure is able to simulate the 
thermal fluctuations of the physical quantities upon ex- 
ploring the phase space with a discrete time evolution 
consisting of subsequent MC steps [thus our time unit is a 
MC Step (MCS)]. The thermal averages become averages 
over the MC evolution after a suitable thermalization of 
the system. We use a thermalization time constituted by 
2000 X n MCS (where n is the total number of the spins 
of our system). Furthermore each point of the Markov 
chain must be generated after some MCS from the previ- 
ous one so that each measure becomes uncorrelated from 
the others. Therefore between two measures we wait a 
time longer than the autocorrelation time. In our specific 
case this waiting time is chosen to be 500 x n MCS. An- 
alyzing the time fluctuations of the observables, we have 
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verified that these choices of thermalization and autocor- 
relation time are optimal to have reliable measures. 

Each configuration can be characterized by the 
instantaneous value of the three magnetizations 
("^(•7r,o)j ™(o,7r)i '^(ji) defining a point in a three dimen- 
sional space. Fig. ^a) is a schematic representation of 
the regions where the points corresponding to instanta- 
neous configurations should be more dense in each one 
of the phases summarized in Table [TTl 

In order to characterize the different phases we com- 
pute the threedimensional (3D) density distribution of 
points at each temperature of the simulation. This dis- 
tribution can be visualized either by projection on one 
plane (Fig.[3Kc)) or by plotting 2D isosurface with a given 
number of points per unit volume (Fig.[3jb)). The points 
are expected to be distributed with a Boltzmann weight 
exp [— F(m(7r^o) > ^'^(o.7r)i thus higher density of 

points corresponds to the lower free energy F and in- 
dicates the most stable phase. As it will be discussed 
in greater detail in the next section, the example in 
Fig. [3Ib,c)orresponds to the nematic phase. Since the 
distribution are rather flat for a given sign of the nematic 
order parameter, the two-dimensional plots like the one 
of Fig. (He) allows to visualize the competition between 
the nematic and the ordered phase. Similar information 
can be cast in the form of a scatter plot as in the exam- 
ple of Fig. [S] below. Notice that the competition between 
the nematic and the disordered phase are not well char- 
acterized by this method. In any case we have checked 
all our results with the more rigorous 3D visualization 
method. Specific details for each phase are discussed in 
next Section. 

Due to the symmetries in the Hamiltonian for a generic 
point (m(7r,o) 1 ''^(o.7r)j ™(/>) there are 7 other points which 
correspond to distinct configurations with the same en- 
ergy. The symmetries are the reflection around the 
'^(ir,o) = s-nd m(o,ir) = planes and the reflection with 
respect to the — plane followed by a 90° rotation 
respect to the axis. This last rotation is due to the 
fact that the sign of is linked to the direction of the 
stripes. With these symmetries for each point obtained 
in the simulations we obtain 7 other symmetry-related 
points increasing the statistic and producing perfectly 
symmetric distributions as expected for non symmetrized 
simulations performed for very long times. 
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Figure 3: (Color online) Phase space of the order parame- 
ters, (a) the spheres or ellipsoids schematically describe the 
different phases of the system. When the order parameter is 
distributed inside the central sphere the system is disordered. 
An order parameter inside the hghtly shaded (orange onhne) 
ellipsoid indicates robust nematic order, while when the order 
parameters are inside the darker shaded (blue online) ellip- 
soids a magnetic state is reaUzed. (b) Example of a MC sim- 
ulation on a 50 X 50 lattice for 7? = 0.51 and fcsT/Ji = 0.52. 
Statistically independent points where binned according to 
the value of the 3D magnetization (m(^,o)j 'Ti(o,7r) , "^0) in cells 
of size 0.1 X 0.1 X 0.1. We show an isosurface with 500 points 
per bin. (c) Contour plot of the distribution projected onto 
the m(o,,r) = plane with bins of size 0.1 x 0.1. The bins with 
higher counts are the most bright. 



netic transition is of first-order—, since the transition is 
'weakly' of the first-order, we can still compute "critical 
indexes" describing the substantial growth of the vari- 
ous susceptibilities near the transition. We first calcu- 
late the temperature and size dependence of the Binder 
Cumulant^ 



U4 [to] = 1 - 



III. RESULTS 

To characterize the transitions to the ordered states, 
we determine the critical indices of the magnetic 
and nematic transitions using the Finite Size Scaling 
techniqueS^i^. Using our MC simulation, we calculate 
the temperature dependence of the interesting physical 
quantities of the system for three different lattice sizes: 
10 X 10, 24 X 24, 50 x 50, using a value of R not too 
close to the T — critical point R = 0.50. Thus we 
choose R = 0.60. Although for this value of R, the mag- 



3(to2 



(5) 



where (m) is the thermal average of the parameter we are 
interesting in (in our case the staggered magnetizations 
or the nematic one) . The temperature at which the value 
of U4 is the same for the three different lattice size defines 
the critical temperature of the transition. We show in 
FigH] the temperature dependence of U4 calculated for 
(to^), and for {mstagg), that is the mean value between 
(to(^ 0)) and (to(o^-)). Our dimensionless temperature T 
is defined in units of |Ji| and taking a unit Boltzmann 
constant fc^ = 1. 
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Figure 4: (Color online) Temperature dependence of U4 [m^] 
and C/4 [mstagg], calculated for different (50 x 50, 24 x 24, 
10 X 10) system sizes. 



From the plots in Fig. 01 one sees that the criti- 
cal temperature of the magnetic transition is Tap — 
0.9715 ± 0.0005, while the critical temperature of the 
nematic transition is T„em ^ 0.9720 ± 0.0005. There- 
fore, within our accuracy, the two transitions occur at 
the same temperature. 

To establish the universality class of the two transitions 
we apply the scaling analysis to the magnetic suscepti- 
bilities 



X(7r,0) 



X{0,7t) 
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N 



N 



(^CT^CTj (-1) 



and to the nematic susceptibility 




(6a) 



(6b) 



(7) 



We consider temperatures sufficiently close to the crit- 
ical temperature, so that we can ignore the second term 
in the r.h.s. of the two equatioirs, {i.e. the square of the 
order parameter) . Indeed at the critical temperature this 
term vanishes like (where t = {T — Tap) /Tap) and 
is negligible with respect to the correlation term (the first 
term in the r.h.s.) that tends to diverge like t~'^ . 

We report the scaling function x = L~^/'^x{t) as a 
function of L^^'^t for the three different system sizes. The 
susceptibilities for systems of different size must coincide 
if the critical exponents 7, v and the critical temperature 
are exact and the transitioir is of the second order. In 
Fig 15] we show the scaling functions for the magnetic (we 
report only X{o.-r)) and nematic suscetibilities. 

We find that a reasonable scaling is obtained with 7 = 
7/6, V = 2/3, Tap = 0.9715. These critical exponents 




Figure 5: (Color online) X(o,7r) and x<t> scaling functions for 
lattice size 10 x 10 (triangles, red and magenta online), 24 x 24 
(circles, green and ligth blue), 50 x 50 (squares, blue and 
black). The temperature step is 5 • 10"*. We have fixed 
7 = 7/6, V = 2/3, Tap ~ 0.9715 for the magnetic transition, 
and 7 = 7.00/4.15, u = 0.975, T„em = 0.9720 for the nematic 
transition. 



correspoird to a four-component Potts model^^. Indeed, 
the magiretic order of our system is described by the two 
parameters TO(7r,o) and TO(o,7r) j each parameter having two 
possible values. Therefore the magnetic order parameter 
has 4 "colors" and threrefore one can naturally expect 
that the transition belongs to the four-component Potts 
model universality class. We remark that the transition 
is weakly first-order-2i and therefore the scaling found can 
only be approximate. 

On the other hand it is natural to expect that the crit- 
ical indexes of the nematic transition are those of the 
simple Ising modeP. We find 7 = 7.00/4.15, v = 0.975, 
Tnem = 0.9720, which are indeed rather close to the crit- 
ical indexes of the twodimensional Ising model, but they 
do not quite coincide with them. Moreover the scaling 
functions of the three different system sizes do not col- 
lapse and are very close only when the reduced tempera- 
ture t is close to zero. Therefore the scaling found in this 
case is not ideal. This is not surprising given that this 
transition is also weakly first-order (see below). More- 
over the two transitions occur simultaneously and can 
influence each other changing the universality class. In a 
more rigorous approach the Ising and the Potts order pa- 
rameters would form a single order parameter with more 
components and a combined universality class. However 
the scaling found above shows that the decomposition in 
Ising and Potts is a good approximation. 

With the present choice of i? = 0.60 we find that there 
is no irematic phase. For example in FiglH] we show the 
scatter plot for the system with size L = 50 at the tem- 
perature T — 0.9720, below which the nematic order 
parameter develops. There is no evidence of a nematic 
phase because there is no accumulation of points in the 
region with q) — and < 0. Notice that the 
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Figure 6: (Color online) Scatter plot in the 2D-spaces 
["^(ir,o) 1 "^0] resulting from the simulations of a 50 x 50 lattice 
with R = 0.60, at the temperature T = 0.9720. 

accumulation of points for > corresponds to (0,7r) 
magnetic order. We find a qualitatively similar behav- 
ior at all intermediate temperatures up to the disordered 
temperature. 

We now consider MC simulations for a 50 x 50 system 
with R = 0.51, that is very close to the critical point 
at i? = 0.50. This case has already been shown as an 
example in FiglH The panel (b) shows a coexistence of 
disorder and nematic order at T = 0.52, because we see 
a maximum in the disordered region and a maximum 
in the nematic region ([?n(7r,o) ~ 0, ^ 0]. Making 
other histograms with higher level surfaces we see that 
the maximum of the disordered phase is lower than the 
one occurring in the nematic region. Therefore the disor- 
dered phase is metastable. The occurrence of two max- 
ima in the phase space population, that is two minima 
in the free energy, also indicates that the transition from 
the disordered phase to the nematic phase is of the first- 
order. Figl3^c) is a contour plot of the counts for each 
bin in the [TO(.n. o),m0] 2D-space. In principle we should 
also study the [m(Q_^), projection. However, since we 
made the two staggered magnetizations equivalent by ex- 
ploiting the symmetries of the problem, we can focus on 
one of the two projections only. 

Decreasing the temperature there is an increase of 
the counts (i.e., configurations) related to the magnetic 
phase, that is identified by finite values of 771(^ 0) ■ The 
temperature at which these counts become higher than 
the ones at nematic region in the 3D representation, cor- 
respond to the critical temperature Taf of the nematic 
to magnetic transition. We show in Fig[7]the level curves 
of the bin counts in the [mT^^o, m^] 2D-space, for a range 
of temperatures. Since also here there is coexistence of 
maxima in the nematic and the magnetic region the tran- 
sition is first-order. 

To compute the temperature range of the nematic 
phase we study the histograms of the simulations at dif- 
ferent temperatures (we use AT — 0.01 steps in tem- 
perature) . A close inspection of the level surfaces of the 
histograms in 3d-space [like the one shown in Fig|31Ib)], 
we can state that at T = 0.53 the maximum in corre- 
spondence to the disordered phase is higher than the 



T = 0.5I T = 0.50 




-1.0 -0.5 O.Q 0.5 1.0 -I.O -0.5 0.0 0.5 1.0 



Figure 7: (Color online) Level curves of the bin counts in 
the 2D-space q) i n^-i/)] resulting from the simulations of a 
50 X 50 lattice, with R — 0.51. Are shown the results for 
T = 0.53, T = 0.52, T = 0.51, T = 0.50. The higher is 
the brightness of the bins and the higher is the values of the 
counts. 




Figure 8: (Color online) Magnetic configurations of a 50 x 50 
systems with T = 0.52; (a) real spin configuration, where the 
red squares represent the spin-up, while the blue squares the 
spin-down, (b) staggered configuration along y axis, where 
the squares are red if at ■ ( — 1)^' ~ +1, while they are blue if 

a, ■ i-iy^ = -1. 

one corresponding to the nematic phase. Similarly at 
T ~ 0.50 the maximum of the smectic phase is the high- 
est. Hence the range of temperatures in which the system 
has nematic order only includes the T = 0.51 and 0.52 
temperatures. We show in Fig|5] the magnetic configu- 
ration at the end of the simulation at the temperature 
T = 0.52. One can clearly see that rotational symmetry 
is broken, but the system does not display a complete 
magnetic order, because the average value of the stag- 
gered magnetization, which define the magnetic order is 
small. We also studied smaller systems with 24 x 24 and 
10 X 10 lattice sites. Performing the same analysis as in 
the previous case, we observe that the range of temper- 
atures where nematic order occurs is larger than in the 
50 X 50 system. We report in FiglS] the same plots as 
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in FiglT] for the 24 x 24 system at the two temperatures 
dehmiting the nematic range. At the temperatures con- 



T = 0.54 
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T = 0.49 




we expect that elongated domains are formed with ap- 
proximately the same aspect ratio. This defines two typi- 
cal length scales and ^|| corresponding to two different 
correlation lengths satisfying 



^11 



(10) 



Approaching the Neel temperature Taf from above, the 
correlation lengths should be described by 



\T-Taf\ 



(11a) 



Figure 9: (Color online) Level curves of the bin counts in 
the [m(.^.o) I "T'.^)] 2D-space resulting from the simulations of a 
24 X 24 lattice, with R = 0.51. The results for T = 0.54 and 
T — 0.49 are reported. The bins with highest counts are the 
most bright. 

sidered in Fig|9] there is still an evident nematic phase. 
For a 10 X 10 lattice, the nematic range is even larger. 
We can therefore conclude that the smaller is the system 
size and the stronger is the tendency to form a nematic 
phase. 

To explain this result we compute the energy cost to 
form a domain like the one displayed in Figl5]on the per- 
fect (tt, 0) order. Notice that this corresponds to a change 
of phase of the magnetic order but without rotation of 
the direction of the stripes. These domains will be nat- 
urally elongated and we aim to find the typical aspect 
ratio. 

We consider for simplicity a rectangular domain like 
the one shown in FiglJJ It is easy to see that the bound- 
ary energy is anisotropic. With respect to the ground 
state (complete smectic phase), a spin located along the 
vertical walls brings along an increase of energy per unit 
length of domain boundary AEy^.^ = 2Ji -I- 4J2- On the 
other hand, the increase of energy per unit length due 
to a spin located along the horizontal wall is given by 
AEhor = — 2 Ji + 4J2. We call l± the length of the walls 
perpendicular to the stripes, and the length of the walls 
parallel to the stripe. The number of spins along the hor- 
izontal walls of the rectangular domain will be 21 ± , while 
along the vertical walls they are 2Z||. Thus the boundary 
energy of the domain is given by 



E, 



bnd 



2?_l(-2Ji +4J2) + 2?||(2Ji +4J2). (8) 



We now determine what is the most stable configuration 
for this type of domains. To this purpose we minimize 
Ebnd varying its aspect ratio while maintaining its 

area l±^l\\ fixed. We find. 



li_ _ 2 J2 + Ji _ 2i? + 1 
77 ~ 2 J2 - Ji ^ 2i? - 1' 



(9) 



Although this computation applies to the ordered phase, 
as the temperature is lowered from the disordered phase. 



el? 



\T-Taf\ 



(lib) 



where and are simple proportionality factors, whose 
ratio is expected to be given by Eqs.® and (|10|) . 

For values of R near 0.5, we have ^ and there 
will be a range of temperatures for which the linear sys- 
tem size, i, satisfies 



Cii < i < 6 



(12) 



This means that the system is ordered in the direction 
perpendicular to the stripes but it is still disordered along 
the direction of the stripes. This is precisely the nematic 
phase. Indeed we can see in the snapshot of Fig. [8] (b) 
that the domains are elongated and span one linear di- 
mension of the sample, while the system has a correlation 
smaller than L in the direction parallel to the stripes. 

As the system size increases the temperature range in 
which Eq. ([T^ is satisfied decreases, i.e. the nematic 
temperature window decreases with system size, in agree- 
ment with our simulations. 

Simple arguments show that the value of R in FeAs 
planes is actually close to its critical value R — 0.5^-. 
These results can also shed light on the experiments by 
Jesche et aL'^ and Liang et ai. The former found that the 
splitting between the structural and the magnetic transi- 
tion decreases as the sample quality is increased while the 
latter finds that the anisotropy in transport decreases. It 
is quite natural to assume that grain boundaries, dislo- 
cations and other defects will disrupt the perfect crystal 
order and introduce an extrinsic scale L roughly given 
by the average distance among defects. As the sample 
is cooled from the disordered phase, at some point the 
longitudinal correlation length will exceed the scale L. 
Then at this scale L the C4 symmetry of the crystal will 
be spontaneously broken and the system will undergo a 
spontaneous orthorhombic distortion possibly with twin- 
ing or domain formation at the scale L. In the perpen- 
dicular direction the system will remain essentially dis- 
ordered preventing the appearance of an elastic signal in 
magnetic neutron scattering. At a lower temperature the 
system will become magnetically ordered in both direc- 
tions. 
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This picture applies for weak magnetoelastic coupling. 
As we will show in the next Section, when the the spin 
system is strongly coupled to the lattice the structural 
and nematic transitions occur at the same temperature 
and the splitting between the structural and the mag- 
netic transition coincides with the splitting between the 
nematic and the magnetic transition. 

IV. SPIN-LATTICE COUPLING 

In this Section we consider lattice deformations and 
their effects on the spin system. In the continuum the 
distortions can be characterized by the symmetric strain 
tensor 

where fi, v label the Cartesian directions. 

In a 2D system it is convenient to work with the fol- 
lowing strains, 

Ci — ^xx ^yyi ^2 — ^xx ^yyi ^3 — V^W^jy, (14^ 

which describe dilational, deviatoric and shear deforma- 
tions respectively. 

To include magnetoelastic effects in our frustrated 
Ising model, we assume that the coupling constant Ju^u 
among atom i and the nearest-neighbour atom in the v 
direction depends linearly on the deformation according 
to 

Jit,x = Ji + "^auxx = Ji + a{ei + 62), (15a) 

Jii.y = Ji+ 2auyy = Ji+ a{ei - 62). (15b) 

For simplicity we will restrict to uniform strains which 
is enough to describe the observed structural transition. 
At lowest order the shear strain does not couple with 
the magnetism. The dilation strain ei simply adds to Ji 
and describes symmetry preserving changes in the vol- 
ume which are interesting but not the focus of the present 
work. Therefore we set ei = 63 = and consider only 
e = 62- The latter is the order parameter of the observed 
structural transition^iii which for a FeAs layers implies a 
deformation from a square lattice to a rectangular lattice. 
Within this approximation we obtain the magnetoelastic 
term 

N N 

Hm-cI = -ae^cri{ai+x - cn+y) ^ -2a£^ </>,;. (16) 

1=1 1=1 

In addition we must consider the elastic energy which is 
given by 

Hei = ^Nfi„e^, (17) 
where fj,o is a shear modulus at T = 0. 



We thus obtain a sum of three terms: the purely mag- 
netic term Eq. ([T]), a purely elastic term Eq. ([T7]) and 
the magnetoelastic term Eq. (|16|) . This last term con- 
tributes to the internal energy U = {H), by an amount 
{HM-ei) — —1Nae{m^) , directly involving the thermal 
average of the nematic magnetization (m^). Thus the 
nematic and structural order parameter couple linearly 
with each other. This implies that the nematic critical 
temperature should coincide with the critical tempera- 
ture of the structural transition. 

To identify the stable equilibrium states of the system 
at finite temperature, we now study the Helmholtz free 
energy. We start from values of temperature higher than 
both the critical temperature Taf of the magnetic tran- 
sition and the critical temperature T5 of the structural 
transition. Therefore no magnetic order and no struc- 
tural deformations are present and the free energy will be 
minimum for (m^) — and e = 0. If we consider small 
variations of the two parameters (m^) and e around their 
zero values, we can expand the free energy up to second 
order in the two variables thereby obtaining 

^^^y^'^ = l{xl)-\m,f + 2a{m,)e+lf,e^ (18) 

Adding external fields coupling linearly with the order 
parameters one can see that is the same nematic sus- 
ceptibility appearing in Eq. ([7]). The last coefficient /i is 
the shear modulus at finite temperature and in the ab- 
sence of magnetoelastic coupling. Within our approxima- 
tions we can consider this term equal to the coefficient fj,Q 
of Eq. (|17p , since the elastic term ([T7| of the Hamiltonian 
coincides with the elastic free energy F^i = fjLe^ . In- 
deed, since the deformation e is constant along the whole 
system, the lattice entropy is a quantity of order 1 and 
can safely be neglected when extensive quantities of or- 
der N are considered. Minimizing Eq. (jl8p with respect 
to the nematic magnetization (m^), we obtain 

{m^)=-2aexl. (19) 
Replacing this expression of {m^) into Eq. (|18p we find 

5F[e) = iiV (mo - 4a\0) = iiY^tot e\ (20) 

where fitot{T) is the temperature dependent shear mod- 
ulus, which takes into account the renormalization due 
to magnetoelastic coupling. 

Introducing an effective coupling constant A — Aa^ / ^q, 
from Eq. (|20p one obtains a generalized Stoner criterium. 
If x^{T) < 1/A, the system is stable, whereas, if x^{T) > 
1/A, the system becomes unstable and gains energy by 
making a deformation s. The temperature at which 
X^{T) — 1/A defines the critical temperature Tj'"* of 
the (second-order) structural transition. 

Let's assume that in the absence of magnetoelastic cou- 
pling, the transitions are second order and that the ne- 
matic and magnetic susceptibility diverge at the same 
temperature, i.e. that there is no true thermodynamic 
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nematic phase for a = 0. The above argument shows 
that for small magnetoelastic coupling the nematic and 
the orthorhombic orders occur simultaneously at the tem- 
perature TJ"'* higher than that of magnetic order. Indeed 
in the disordered phase, the nematic susceptibility 
will grow upon lowering the temperature and the insta- 
bility (Stoner-like) condition will be satisfied while the 
magnetic susceptibility is still finite. Although the pres- 
ence of the structural distortion will tend to increase the 
Neel temperature T^p, this increased temperature can 
never reach Tj"'' since e vanishes at Tj"'^. Therefore a 
finite temperature window should appear between Tj"'* 
and Taf- Thus the magnetoelastic coupling can induce 
a nematic phase. Of course, if the nematic phase were al- 
ready existing for a = 0, a small magnetoelastic coupling 
would tend to widen the nematic temperature window. 

Notice that the above arguments are general and do 
not depend on the details of the model (e.g. Ising vs. 
Heisenberg) but on the second order character of the 
transition. We will show below that the weakly first-order 
character is enough to change this scenario. For a first 
order transition occurring at T^^^ the temperature Tj'"* 
where the shear modulus vanishes satisfy rj'"^ < Tg^^ 
and corresponds to the limit of stability of the metastable 
disordered phase. 

In order to investigate the feedback effect of the lattice 
on the magnetic transition and the effect of the first- 
order character we carry out a Monte Carlo study of the 
magnetoelastically coupled Hamiltonian. To take into 
account the thermal fluctuation of the deformation pa- 
rameter AJ = ae, every n Monte Carlo step (where n 
is the number of the lattice sites) we allow for random 
variations of this parameter in a certain range (we use 
[—0.2 : 0.2]). Each variation of A J is accepted or re- 
fused according to Metropolis algorithm with the energy 
given by Eqs.(IT|), and ([T7]). We perform temperature 
runs (with a step AT = 0.01) of a 50 x 50 system with 
R = 0.51, for several values of A. 

The phase diagram in the [A, T] space, resulting from 
our simulations, is reported in FiglTUl We see that in- 
creasing A, the temperature interval with nematic order 
shrinks and it disappears completely for A — 2.0 • 10^^. 
Thus the largest nematic window is for A = and it cor- 
respond to the case studied in the previous section. This 
goes against our original expectation to find an enhanced 
nematic phase in the presence of the magnetoelastic cou- 
pling. Such a failure is clearly due to the first order char- 
acter of the transition which gets enhanced by the mag- 
netoelastic coupling. This can be recognized from the 
fact that rj"'* and T^"* practically coincide for A = 0, 
while Tg''* becomes substantially larger than Tj"'^ as A 
increases. This is also supported by the presence of ther- 
mal hysteresis of the order parameters that we observe 
at high values of A. 

Notice that as expected Tj"'* increases with A but Taf 
increases more and due to the first-order character of 
the structural transition this renormalization of Taf can 
reach Tj^* closing the nematic window. As discussed 
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Figure 10: (Color online) Phase diagram of the system with 
R = 0.51 in the A vs. T] space. We report the numerical 
results for the structural (or nematic) critical temperatures 
(open circles, blue online) and the magnetic critical tempera- 
tures (open squares, green online), and the structural critical 
temperatures in the stability limit (open circles, red online, 
joined by the solid line). The shaded (magenta online) region 
represents the nematic range in a 50 x 50 lattice. The inset 
displays an expansion of the nematic region of the phase dia- 
gram. The nematic range in a 24 x 24 system (darker shaded 
region, yellow online) is also shown. 



in the previous Section, this window is larger for smaller 
system sizes. Indeed, as it is shown in the inset of Fig lTUl) . 
in a 24 X 24 system a larger nematic region in the [A,T] 
space is found. Although increasing A leads to the closing 
of the nematic window, a nematic phase does exist at 
small spin-lattice coupling: We report in Fig[TT]a contour 
plot of the order parameter distribution function for A = 
1.0 • 10~^ which shows clearly that the system is in a 
nematic phase. 




-0.5 0.0 0.5 
'"(1,01 

Figure 11: (Color online) Contour plot of the order parameter 
distribution projected in the [m(^_o)i ^4>] plane resulting from 
the simulations of a 50 x 50 lattice, with R = 0.51 and A = 
1.0 ■ 10"^ at T = 0.538. 

In short we can state that for small A, transitions 
from disordered phase to nematic phase and from ne- 
matic phase to magnetic phase are both only 'weakly' 
of the first-order. Therefore our schematic model sys- 
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tern grossly reproduces the features of the 1111-3^ or 
of Ba(Fei_a;COa;)2As2^iiS, materials, where a separa- 
tion between a second-order structural transition and 
a second-order magnetic transition is observed. In 
particular we notice that, upon increasing doping in 
Ba(Fei_a;Coa;)2As2 the nematic window increases becom- 
ing largest at a; > 0.06, near optimal dopingl^. This 
strengthening of the nematic state finds a rational within 
our findings because increasing doping naturally leads to 
a decrease of magnetoelastic couping A and to an increase 
of disorder. As found in our schematic model, both these 
effects increase the nematic region between the structural 
and the magnetic transitions (cf. the inset of Fig. [TO)) . 
We also notice that estimates from shear modulus mea- 
surements in Ba(Fei_a;Coa;)2As2 give typical values of 
A < 5 • 10"'^, substantially smaller than the values at 
which the nematic windows closes in our Fig. 1101) . 

When instead we choose a high value of A, we obtain 
a system which grossly reproduces the behavior of many 
122 materialsii where the structural and magnetic transi- 
tions are strongly first-order and occur at the same tem- 
perature. We furthermore notice that when the tran- 
sitions become strongly of the first-order, the value of 
the critical temperature increases. Although the precise 
value of the critical temperature will depend on many 
other parameters it is interesting that the critical tem- 
perature of the '122' materials are substantially higher 
than the critical temperatures of the '1111' materials. 



V. CONCLUSIONS 

In this paper we investigated numerically and with 
simple analytic arguments the possible occurrence of a 
nematic phase and the associated splitting of the struc- 
tural and magnetic transitions in FeAs planes. We ana- 
lyzed a schematic model of Ising spins with nearest and 
next-nearest neighbor couplings. The effects of a spin- 
lattice coupling has also been considered. First of all 
we find that for R larger than (but close to) 1/2 the 
model displays a weakly first-order transition. The very 
weak first-order character found in the absence of lattice 
coupling allows to study critical exponents which for the 
magnetic order parameter, are close to the ones of a 4- 
state Potts model in agreement with simple symmetry 
considerations. 

One can then remark that the magnetic Potts order 
parameter can order at finite temperature even in two 
dimensions and therefore long range magnetic order ef- 
fectively competes with the Ising nematic order. This 
model, therefore, has the double advantage of being more 
easily solved in numerical MC approaches and of being 
an "acid" test for the occurrence of nematic order. In 
addition, as we argue in the introduction, collective mag- 
netism produces a system with large connectivity of the 
magnetic Hamiltonian, which makes the thermal transi- 
tion rather mean-field like. We argue that such mean- 
field transition in a real lattice with 3D couplings and 



weak anisotropics is more similar to a 2D Ising transi- 
tion than to a 2D Heisenberg transition. Thus, despite 
the apparent obvious diff'erence in symmetry between our 
model and a real system, the model for some aspects is 
expected to be more realistic than a frustrated 2D Heisen- 
berg model. 

Our numerical study indicates that in the region of 
the stripe ground state, but close to the transition to the 
Neel AF phase {R > 0.5), nematic states can form due 
to finite size effects in some temperature range above the 
magnetic order (the smaller is the length scale and the 
larger is the temperature range). We argue that such 
finite size effects are a proxy for the effect of crystal im- 
perfections which disrupt the perfect lattice order and 
provide a natural explanation for the results of Jesche 
et al* and Liang et al.^ who find more robust nematic 
signatures in the worst quality samples. 

The increase in resistivity anisotropy^ with decreasing 
sample quality has also been explained by Fernandes et 
ai— in terms of the interplay between scattering of car- 
riers by spin fluctuations and disorder.^ We believe both 
Ref. i9| mechanism and the present one can contribute to 
the observed effect. More experimental and theoretical 
work would be needed to determine which one will be 
dominant in the different regimes. 

We have neglected all Heisenberg physics which is ex- 
pected to enhance the tendency to nematic phases i^ii^ 
However our arguments for the enhancement of the ne- 
matic phase in imperfect systems also applies to other 
models. For example also in a more realistic 3D Heisen- 
berg model, possibly with long range couplings and 
anisotropies, the domains will be elongated and there will 
be two characteristic length scales, so that our arguments 
remain valid. 

It is sometimes argued that, since the structural tran- 
sition occurs at a higher temperature than the magnetic 
transition, it is the structure that drives the magnetism. 
This argument is clearly wrong and our model is another 
counterexample which show that the structural transition 
can occur at a higher temperature than the magnetism 
although it is driven by the magnetism itself. This is be- 
cause the structural transition is driven by the cohesive 
energy which is determined by short range correlations 
like {(TiCTj) with i and j close neighbors. Thus cooling the 
system one can have a robust short range (aiaj), which 
renders the lattice unstable, before the system has long 
range magnetic order. Formalizing this argument in the 
case of second order transitions, one obtains the addi- 
tional result that the coupling to the lattice favors the 
nematic state. However in the simulations we find the 
opposite because the lattice tends to make the transi- 
tion more first-order like. This suggest that the nematic 
state observed is mainly driven by the electronic degrees 
of freedom and that the lattice act as a spectator, as 
suggested in Ref. i6j. 

The different behavior found for weak and strong mag- 
netoelastic coupling corresponds well with the behavior 
found on different pnictide families. Indeed 122 mate- 
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rials have larger critical temperature, first-order transi- 
tions a non-nematic state, as found for strong magnetoe- 
lastic coupling. On the other hand 1111 systems have 
lower transition temperatures, second-order or weakly 
first-order transitions and a nematic state as found for 
weak magnetoelastic coupling. 
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